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Abstract. A theory of stress fields in two-dimensional granular materials based on directed force chain 
networks is presented. A general Boltzmann equation for the densities of force chains in different directions 
is proposed and a complete solution is obtained for a special case in which chains lie along a discrete 
set of directions. The analysis and results demonstrate the necessity of including nonlinear terms in the 
Boltzmann equation. A line of nontrivial fixed point solutions is shown to govern the properties of large 
systems. In the vicinity of a generic fixed point, the response to a localized load shows a crossover from 
a single, centered peak at intermediate depths to two propagating peaks at large depths that broaden 
diffusively. 

PACS. 45.70.Cc Static sandpiles; granular compaction - 83.80.Fg Granular solids 



1 Introduction 



The response to a localized force applied at the bound- 
ary of a semi-infinite sample is an essential feature of 
any macroscopic material. For materials that are well de- 
scribed by linear elasticity theory, the response can be 
calculated by standard techniques: using the relation of 
the stresses to a displacement field, one constructs partial 
differential equations for components of the stress tensor 
cr, then solves them with appropriate boundary conditions 
on a and its derivatives. For a wide class of granular ma- 
terials, however, the absence of an energy expressible in 
terms of microscopic (or grain scale) displacements leads 
to serious difficulties in deriving the stress response. 

Numerous attempts have been made, many quite re- 
cent, to close the system of stress equilibrium equations 
by deriving (or simply guessing) a relation between the 
components of the stress tensor beyond those required 
by Newton's laws. The needed relation may be extracted 
from assumptions concerning yield thresholds in elasto- 
plastic theories [|l|j2|, from an analysis of physics at the 
grain scale, |^ , from general symmetry principles and con- 
siderations of simplicity , or from considerations ap- 
plicable to isostatic networks A special case that 
has received much attention is that of frictionless, circular 
disks. 1^,^ Alternatively, lattice models for the configura- 
tion of individual inter-grain forces can serve as the basis 
for the derivation of average stress patterns while simul- 
taneously addressing the statistical properties of forces at 
the grain scale. 



In this paper we develop a theory of stress distribution 
in noncohesive granular materials based on the physics 
of force chains, rather than macroscopic stresses or grain 
scale interactions. Experimental probes and numerical com- 
putations of stress patterns in granular materials univer- 
sally show filamentary networks that support strong forces 
between grains. [ p^l|l^ , [l5| We take this as motivation to 
describe all of the forces in the system in terms of linear 
"force chains" that combine to form a network whose large 
scale properties determine the macroscopic stress field. To 
describe such networks, Bouchaud et al. recently intro- 
duced the "AY- model" (pronounced "double Y model" ) . 
Here we develop the AY-model in greater detail, point out 
certain essential ways in which the original formulation 
and analysis was incomplete, and calculate stress patterns 
and response functions in special cases. 



In the AY-model, stresses are modeled as networks of 
interacting line segments that carry compressive forces. A 
crucial feature of the model is that the force chains are 
directed: each is conceived as "propagating" from one of 
its endpoints to the other, in a sense that will be clarified 
below. Chains are created by three processes: (1) imposed 
conditions at a boundary; (2) the splitting of one chain 
into Z) in a manner that preserves force balance at the 
splitting point, where D is the spatial dimension of the 
material; or (3) the fusion of D chains into one. They 
then propagate until they are destroyed via splitting ( "A" 
processes) , fusion ( "Y" processes) , or traversal of a sam- 
ple boundary. We refer to a static network of chains in- 
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teracting in this way as a "directed force chain network" 
(DFCN). 

Our goal is to determine the large scale structure of 
DFCN's and their average responses to small, local per- 
turbations at a boundary. We first develop a Boltzmann 
equation governing the spatial variations in the densities 
of chains supporting certain force intensities and oriented 
in particular directions]^ We then seek solutions in the 
special case where the chains are restricted to lie along a 
discrete set of directions. The results highlight several sub- 
tle features of the directed force chain system, and show 
that the response of such systems can have a rich structure 
that includes some surprising eff'ects. 

Our primary result can best be stated with reference 
to the results previously obtained via standard elasticity 
theory or posited closure relations. In standard elastic- 
ity theories of two-dimensional isotropic systems, the re- 
sponse to a normal force applied at one point (x = 0) of 



a half plane is 1 1 7 
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where z = at the surface and increases with depth. 
Thus at any fixed depth z, the response CTzz consists of 
a single peak centered on a; = 0, with a width that grows 
linearly with depth. Though a^z for anisotropic materi- 
als can have two peaks, the response in this geometry is 
always invariant under uniform dilations and rescaling: 
(x ax, z az) yields try {l/a)aij. 

In another class of theories based on closure relations 
that make no reference to a displacement field, the stress 
equilibrium equations are hyperbolic. In such cases, the 
response function consists of two peaks that propagate 
linearly away from a; = 0. The inclusion of weak disor- 
der in these models leads to a diffusive broadening of the 
peaks with increasing depth. ||^ Note that this form of re- 
sponse differs markedly from the elastic the ra- 
tio of peak width to propagation distance decreases like 
1 / -^/z at large depths rather than remaining constant. Per- 
haps even more importantly, the two types of theory re- 
quire different manners of specifying boundary conditions. 
Whereas standard elasticity theory requires the specifica- 
tion of two conditions on the stresses (and/or their deriva- 
tives) at all points on the boundary, the hyperbolic models 
require (permit) the specification of stresses on only the 
top boundary. 

In the directed force chains system, we find that both 
top and bottom boundary conditions are important, but 
the solutions and response functions in deep systems have 
a richer structure than either the elastic or hyperbolic the- 
ories suggest. For very shallow depths beneath the surface, 
the response may have one or two peaks, depending on the 
details of the applied forces. For intermediate depths, up 
to several times the average length of a force chain, the 
response has a single peak and may appear quite similar 
to a standard elastic response. For large depths, however. 



^ We use the term "Boltzmann equation" loosely. It does not 
imply the existence of a causal structure or H-theorem. 



two peaks emerge that grow diffusively with depth. This 
latter behavior suggests that at the largest scales a hy- 
perbolic model for stress equilibrium may be appropriate 
even for systems with strong disorder, though important 
questions remain concerning the role of the discrete set of 
directions in sustaining the two-peaked structure. 

A second result to be emphasized is the failure of the 
linearized theory obtained by assuming that all chain den- 
sities are small. Formally, this assumption is equivalent to 
neglecting the fusion of chains, keeping only those terms 
in the Boltzmann equation that describe splitting events. 
For intermediate and deep systems, internal consistency of 
the theory requires that the response be computed by lin- 
earizing around a nontrivial fixed point of the Boltzmann 
equation. This sheds new light on previous efforts to derive 
a linear elasticity theory of DFCN's. It also has important 
implications for the numerical modeling of DFCN's and 
the interpretation of experimental measurements of stress 
response in granular systems. 

In addition to these two central results, we point out 
several curious features of DFCN's in a simple slab geom- 
etry in two dimensions. These include an effective Pois- 
son ratio that may depend on the way in which a load 
is applied, an exponentially localized region of increased 
horizontal stresses in the middle of a slab subject to hori- 
zontally uniform vertical stresses at top and bottom, and a 
distinction between responses to identical applied macro- 
scopic stresses with different distributions of force chain 
densities in various directions. 

The paper is organized as follows. In the Section ||, 
we introduce the concept of the directed force chain, de- 
fine the constants, variables and functions that enter our 
theory, develop a Boltzmann equation that governs spa- 
tial variations of the densities of force chains, and discuss 
the boundary conditions on the equation. We then define 
a special case of the model in which chain directions are 
restricted to discrete set that is amenable to analytical 
treatment. In Section ^, we solve the discrete Boltzmann 
equation in a slab of depth z for the case of applied loads 
that are uniform in the horizontal direction, identifying 
the (nontrivial) solutions relevant for the computation of 
response functions. In Section |^, we analyze the response 
(in the discrete model) in the vicinity of a generic fixed 
point, present some results for other choices of discrete 
directions in 2D and describe a possible generalization to 
3D. We also discuss the relation of our analysis methods 
to the long- wavelength theory presented in Ref. We 
conclude, in Section o with discussions of the issue of nu- 
merical generation of DFCN's and some open questions. 



2 Boltzmann equation for directed force 
chains 

2.1 Definitions, the Boltzmann equation, and its 
boundary conditions 

A force chain is defined as a line segment that supports 
a compressive force. The direction of a force chain has 
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Fig. 1. Illustration of force chains, (a) Two force chains initi- 
ated at the top surface cross without interacting, (b) When one 
of the force chains is initiated at the bottom, a fusion becomes 
possible, (c) A force chain initiated at the top that splits upon 
encountering a defect. 



two features. First, the segment makes a particular angle 
with the vertical. Second, the two ends of the segment 
are distinguished - the segment has a "beginning" and an 
"end" determined by the role the chain plays in the net- 
work. Fig. illustrates the latter point. The figure shows 
a two-dimensional system consisting of a hexagonal array 
of grains. In panel (a), two forces are applied at the top 
boundary and the forces propagate along chains, crossing 
at the central grain. In panel (b), a similar situation is 
shown, but this time one of the chains is assumed to be 
specified by fixing its position at the bottom boundary in- 
stead of the top. In this fusion of the two chains at 
the central grain is permitted (though not required) and 
the resulting configuration may be different, demonstrat- 
ing the inequivalence of the two possible choices of direc- 
tion for the chain whose boundary condition was changed. 
The horizontal chain is assumed to "begin" at the fusion 
point and "end" at the boundary. 

Thus each chain has an intensity, /, and a direction n. 
Though we will sometimes combine the two and express 
them as a vector f , it is important to keep in mind that 
the compressive force supported by the chain has no for- 
ward or backward direction. A grain that is part of a force 
chain experiences no net force, only a stress with its major 
principal axis along the chain. The distinction between the 
forward and backward directions of the chain refers to its 
own boundary conditions, which are ultimately related to 
the boundary conditions at the edge of the sample, though 
this relation may be difficult to determine. 

In addition to fusion of chains, it is possible for a chain 
to split into two (in two dimensions). Fig. [^(c) shows a 
splitting event induced by the presence of a defect in the 
system. In general, ordered regions and defects will not 
be so easy to identify by eye, and all sites are assumed to 
produce splittings and permit fusions with probabilities 
drawn from a specified angular distribution, as described 
below. 

To clarify the way in which directions may be assigned 
to force chains, it may be helpful to consider the config- 
uration shown in Fig. |^. Here the question arises as to 
how one can assign a direction to the chain in the mid- 
dle of panel (a). If the situation is as shown in (b), with 
boundary imposed at the two sites on the top surface, then 
the middle chain must be interpreted as consisting of two 
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Fig. 2. (a) A seemingly ambiguous configuration of force 
chains, (b) Boundary conditions that imply an annihilation 
event, (c) Boundary conditions that give a unique direction to 
the middle chain. 



oppositely directed chains that annihilate in the middle. 
The same is true for any other configuration of bound- 
ary conditions that force this configuration on the mid- 
dle chain. Alternatively, there exist other possible speci- 
fications of the boundary conditions that give the middle 
chain a specific orientation, as shown in (c). In develop- 
ing the theories below, we neglect all annihilations of the 
type shown in (b). The reason is that for very narrow force 
chains, the probability of two chains meeting head on is 
very low. In granular systems, it may be argued that the 
force chain widths are effectively of the order of a grain 
size, so that annihilation events are not entirely negligi- 
ble. For the present paper, we neglect effects associated 
with the finite sizes of grains. Preliminary investigations 
strongly suggest that inclusion of terms corresponding to 
annihilations does not change the general features of the 
results. 

For clarity of presentation, from here on we work in two 
dimensions unless explicitly noted otherwise. The main 
difference in D dimensions is that vertices where chains 
fuse or split generically require D + 1 chains in order 
to achieve force balance, since force balance among fewer 
chains would require the highly improbable event that one 
of them lie in the subspace spanned by the others. 

We begin with the assumption that a 2D granular ma- 
terial can be thought of as a collection of local environ- 
ments through which single force chains are either trans- 
mitted or split into two, and pairs of intersecting force 
chains either pass through each other or fuse into one. We 
also assume that the force chains carry all of the stresses 
in the system. Newton's laws are built into the model via 
the constraint that the forces at each splitting and fusion 
event must balance. In this formulation of the model, there 
are never any unbalanced torques, as the three chains as- 
sociated with a splitting or fusion always meet at a single 
point. 

Following the suggestion of Ref. we construct a 
Boltzmann equation governing the densities of chains as 
follows. Let P(/, n,r) represent the probability of find- 
ing a force chain of intensity / in direction of unit vector 
fi passing through the spatial point r. In other words, 
J^^^df J dn |n-u|P(/, n, r) is the number of chains of in- 
tensity between g and g + S that cross a unit length line 
segment through r and perpendicular to u. Note that P 
has dimensions of 1/ (force x length) in 2D and cannot 
be negative. The local stress tensor, expressed in terms of 
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the coarse grained quantity P, is given by |If 



/■oc /• 

J df Jdnn^npfP{f,n,r) 



(2) 



Here f2 represents the D-dimensional angular direction of 
n divided by the full solid angle so that / df2 = 1. For our 
2D discussion, we will use 6 € (— tt, tt] to designate the 
direction fi, and d/2 = d9/2TT. 

We define two functions 4's and (f)f that characterize 
the probability of generating chains at various angles when 
a splitting or fusion occurs. In the most general case, these 
may be functions of the intensities and directions of all 
three chains, but must include delta functions ensuring 
force balance and Heaviside functions 0(f) ensuring that 
all force intensities are positive (all force chains in a nonco- 
hesive material must carry compressive stress, not tensile). 

M^i I f2, fa) = sih + fs - fi) e(/i) eih) 0{h) 

>^M0i\e2,0^)\sm{e2- 0:^)1 (3) 

where f indicates the pair (/, n) and tpa is normalized 
such that jd62d6^df2df^4)a — 1. The vertical bar is used 
in expressing the arguments of (pa and ipa to indicate that 
the first argument is associated with the minus sign in 
the delta function. For splittings, the first argument cor- 
responds to the single incoming chain. For fusions, it cor- 
responds to the single outgoing chain. 

The exphcit factor of | sin (6*2 — 6*3)! in Eq. (§) serves 
several purposes. For the splitting function, it is included 
to remind us that 4>s must vanish for splitting events in 
which the outgoing chains are coUinear and oppositely di- 
rected, as this would generate infinite force intensities. It 
is also reasonable to assume that splitting events with out- 
going chains close to the same direction are rare, though 
•ips could of course be chosen to make such events as prob- 
able as desired. For the fusion function, the factor gives 
the density of intersection points per unit area for unit 
chain densities of chains with orientation 62 and 63, which 
should clearly affect the probability of fusions occuring. 
This factor also simplifies the normalization of tpa, since 
for 01 = and any function A we have 

1*00 

/ df2dh5{i2 + h-h)A{hj2,h) (4) 
Jo 

1*00 

= / #3 5 (/2 cos 62 + /3 cos 6*3 - /i ) 

X (5(/2sin6'2-l-/3sin6l3)^(/i,/2,/3) 
1 _/isin6'2 /i sin 6*3 



sin(6'2 - 6*3) V sin(6l2 - 9^) ' sin(6i2 - 6*3) 



Note that ipa must be symmetric in its second and 
third arguments; in the case of splitting (fusion), switch- 
ing the labels of the two outgoing (incoming) chains can- 
not change the physics. Note also that the delta function 
is really the product of two delta functions, as shown in 
Eq. (^, one for each component of force balance, and 
therefore has a dimension of inverse force squared. We 
will restrict our attention to homogeneous and isotropic 



systems, for which ipa does not vary with position and 
depends only upon the relative angles between its three 
arguments. 

In general, ips and ipf are different, owing to the phys- 
ical difference between splitting and fusion events. As il- 
lustrated in Fig. |l| for a particularly simple system, the 
type of defect required to induce splitting is not necessar- 
ily present at a fusion. The boundary conditions on the 
various chains involved make for different probabilities be- 
tween splittings and fusions for a given local geometry of 
the granular packing. 

The following Boltzmann equation describes the vari- 
ation of P{f: 9, r) along the direction n, where A and /i 
are constants required on dimensional grounds and, for 
notational convenience, we have dropped the r from the 
argument of all of the P's: 



(n . V) P(f) = 



P(f) + 2 df'df"d9'd9" (PsiJ' I f, f")^^(f') 



dfdf"d9'd9" [0/(f I f, f")P(f')^(f") 

^24>f{i'\fX)p{wn]- (5) 



In this expression, all integrals over angles run from — tt to 
TT and all integrals over forces run from to 00. The first 
term on the right hand side in Eq. (|^) describes the decay 
of force chain density in the direction of the chain due to 
splitting of the chain. A thus corresponds to the average 
distance a chain propagates before splitting, in the absence 
of all other chains. In principle, A may be a function of 
the direction n and position, but is a single constant for 
homogeneous and isotropic systems. The second term on 
the right hand side describes the increase in P(f) due to 
the splitting of other chains. The factor of two in this term 
is due to the fact that 4>s has been normalized to unity 
while a splitting eventproduces two outgoing chains. The 
final integral in Eq. (^ is quadratic in P and describes 
the effects of fusions. The quantity ^ has dimensions of 

or force x length in 2D. As with A, we assume 
is independent of 9. The factor of 2 in arises because the 
chain in the 9 direction can be either of the two incoming 
chains in a fusion. 

The boundary conditions to be applied to Eq. (^ are 
suggested by Fig. |l]. At each point on the boundary of 
a sample we may specify the density of ingoing chains 
of any intensity, but we may not specify the density of 
outgoing chains. This means that there is no simple way 
to assign boundary conditions corresponding to the over- 
all stress on a boundary; the density of outgoing chains, 
which add to the normal and shear stresses just as in- 
going ones do, is determined by the splittings and fusions 
within the sample. In the slab geometry, however, it is pos- 
sible to specify the total vertical force apphed to the top 
and bottom boundaries, due to the symmetry of the force 
network under multiplication of all intensities by a com- 
mon factor. One may first specify the density of inward 
directed chains, then compute the density of outward di- 
rected chains using the Boltzmann equation, compute the 
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stresses at the boundary, then rescale ah of the forces in 
the system to arrive at the desired boundary stress. Sat- 
isfying conditions on both the normal and shear stress 
may also require adjustment of the relative densities of 
chains in different inward-pointing directions. The man- 
ner in which boundary conditions on the chain densities 
can be specified will be clarified by the precise treatment 
of a special case in Section |[ 

The fact that we may only specify the densities of in- 
going chains at a boundary is a part of the definition of 
the AY-model that deserves further comment. Microscop- 
ically, when we specify that a chain of a given strength 
must pass through a certain point on the boundary, the 
boundary condition defines the direction of that chain to 
be inward. This is physically plausible because we would 
not expect to be able to apply a boundary condition that 
requires an intricate conspiracy of fusions and splittings 
to create a chain that propagates outward at the point 
in question. In the context of the Boltzmann equation, in 
which we do not specify individual chains but only densi- 
ties, this is reflected mathematically in the fact that spec- 
ification of outgoing chain densities leads to nonsense in 
certain generic situations, as shown in Section ^ for ex- 
ample. 

To appreciate the meaning of the chain directions at 
the boundaries, it may be helpful to consider the appli- 
cation of force in the following way. Imagine a packing of 
grains between to flat plates that are perforated with holes 
much smaller than the grain size and held at fixed separa- 
tion. Now imagine poking needles through the holes in the 
plates, with the force applied to each needle being speci- 
fied externally and an equal total force applied from above 
and below. In this situation, there is a clear difference be- 
tween chains that end on needles and chains that end on a 
contact between a grain and one of the plates. The former 
must be present due to the boundary conditions, and are 
therefore ingoing. (See Figure |[) The latter exist only as 
a response to the applied forces. They would shift around 
if a new needle were poked in initiating a new chain and 
hence are not specifiable as boundary conditions . Now in 
the case of two rigid, smooth flat plates that are simply 
pushed together, it is not clear how to distinguish the in- 
going from the outgoing chains. Nevertheless, for present 
purposes, we assume (plausibly, we think) that such a dis- 
tinction is somehow embedded in whatever physics at the 
grain scale is responsible for the organization of stress into 
force chains in the first place. Exploration of the alterna- 
tive hypothesis ~ that ingoing and outgoing chains are 
indistinguishable - is beyond the scope of this paper. It 
would require self-consistent choices of the splitting and 
fusion functions so as to yield structures that satisfy the 
same Boltzmann equation when the directions of all force 
chains are reversed. 



2.2 Rescalings and the importance of nonlinear terms 



ing of lengths: x x/X, z ^ zjX^P ^ XP and /i — > /i/A. 
In these new units, fi has dimensions of force. 

A crucial feature of Eq. (H) is that the fixed point so- 
lution P(f) = is unstable. Consider the total force den- 
sity in the system, defined a,s J- = JdfdddrfP{{). For 
P(f) << 1, the value of T is determined by the linear 
terms in the Boltzmann equation. But these describe only 
splittings, and every splitting increases since force bal- 
ance requires that the sum of the two outgoing intensities 
be larger than the incoming intensity. If we begin with 
a single chain, each generation of splittings increases 
by a constant factor that depends upon ips but is always 
greater than unity. In a sufficiently large sample, the rate 
of increase of J- due to splitting will exceed the rate of 
decrease due to leakage through the boundary. In order to 
regulate this divergence, the nonlinear terms must be in- 
cluded. This argument wil l be made precise for the choice 
of Ips discussed in Section ^.3| and developed in somewhat 
more general terms in the appendix. 

The inclusion of the nonlinear terms requires the intro- 
duction of the dimensionful parameter fi. To understand 
how this parameter is determined, it is helpful to write it 
as the product of two factors, fi = fJ-oY, where /io has units 
of force and Y is simply a numerical constant. Since there 
is no intrinsic force scale in the system, /zq can be chosen 
arbitrarily without loss of generality. This freedom is a di- 
rect consequence of the fact that, given any DFCN, there 
exists a continuous family related to it by multiplication 
of all of the chain intensities by a constant factor. The 
rescaling associated with the choice of /xq is exactly com- 
pensated by a rescaling of the argument / in P{f, 9, r), so 
that the geometric structure of the DFCN is unchanged. 

The factor Y is determined by the overall probabil- 
ity of fusions occurring when chains intersect. Y can be 
thought of as a normalization associated with In phys- 
ical terms, changes in Y have real effects on the P's, as Y 
ultimately determines the density of fusions in the sys- 
tem. It is convenient, however, to work with variables 
P' = P/Y. This results in the same Boltzmann equation 
with Y = 1, which we will work with henceforth while 
dropping the primes on the P's. Thus we are left with a 
Boltzmann equation with no parameters other than the 
choices of ips and ipf. for the purposes of mathematical 
analysis A and /i can both be scaled to unity. 

Eq. (^), with A and fi both scaled to 1, but with generic 
forms of i/js and tpf, has proven difficult to solve for even 
the simplest sample geometries and boundary conditions 
(though the authors still hold onto hope for further analy- 
sis). Moreover, the numerical generation of directed force 
chain networks encounters serious difficulties for systems 
deep enough that the nonlinear processes become impor- 
tant. Useful insights can be obtained, however, by con- 
sidering some special cases where analytical progress is 
possible. 



In the isotropic Boltzmann equation, the parameter A can 
be set equal to unity without loss of generality via a rescal- 
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Fig. 3. Two choices of the discrete set of six permitted direc- 
tions of force chains. 



2.3 A solvable case: 120° splittings and fusions 



We now consider a system for which 



72, C3 



) = 



1 



5(921 - 
+5(021 



l)S(03i + |) 



(6) 



where 9ij = 6i — 9j. This corresponds to the case where 
both sphttings and fusions are always symmetric and form 
vertices composed of 120° angles, as in Fig. |ll Note that 
in this case, chains with different intensities or with direc- 
tions that are not parts of the single symmetric star of six 
unit vectors, cannot interact via fusion. So we further as- 
sume that all chains initiated at the boundary have direc- 
tions chosen from a single "6-fold star" and have identical 
intensities /*: 



(7) 



n=0 



where 9n = (n — l/2)7r/3 as illustrated in Fig. ||(a). We 
refer to this orientation of the 6-fold star vectors as "hor- 
izontal" , since it includes directions oriented horizontally. 
Another choice we have studied extensively is the "ver- 
tical" orientation shown in Fig. ^(b). It turns out that 
for the horizontal orientation, exact analytical expressions 
are easier to obtain, but the response functions (computed 
numerically for the vertical orientation) are qualitatively 
similar. 

For the horizontal orientation, Eq. (j^) reduces to 

V, • VP, = -P, + P,+i + P,_i 

-f (P,+ lP,-l - P,P,+2 - P,P,-2) , (8) 

where all indices are taken modulo 6 and = (cos 9i , sin 9i ) . 
Note that all integrations over / and 9 have been per- 
formed, and the only remaining difference between the 
coefficients of splitting and fusion events is the factor Y, 
which has been absorbed into the P's. 

The restriction of the possible orientations of force 
chains to six directions does not imply that the chains 
reside on a regular geometric lattice. The splitting and fu- 
sion events occur at arbitrary positions, governed only by 
the probabilistic rules encoded in the Boltzmann equation. 



3 Horizontally uniform solutions 

Let us now construct solutions of Eq. (^) for the case of an 
infinite horizontal slab subject to a horizontally uniform 
load. Let z — represent the top surface of the slab and 
z = d the bottom. For simplicity, we consider the case 
where the loading (and hence the full solution) is sym- 
metric under reflection through a vertical line, as well as 
translationally invariant in the horizontal direction. Thus 
we have Pq = Pi, P5 = P2, and P4 = P3 at all points in 
the system and dxPn =0 for all n. 

The Boltzmann equation then reduces to the following 

set: 



V3 



9.P1 



P1P3 



(9) 



P2 - Pi + P3 + [P1P3 - P2(Pl + P3)] , (10) 



V3 



a.Pj - P2 - P1P3 



(11) 



where each of the P's is a function of z alone. From Eqs. (^ 
one sees immediately that G = Pi -I- P3 is a con- 
ctermined by the boundary conditions), and from 
wc have 




P2 = 



G + P1P3 
1 + G 



(12) 



Substituting into Eqs. (p|) and (^TJ) leads to the two cou- 
pled ordinary differential equations for Pi(z) and P3(z): 



dPi 
dz 

dz 



2G 



(1 



-G)V3 
-2G 



(l + G)\/3 



(1 - PiPz) 



(1 - P1P3) ■ 



(13) 



These equations are to be supplemented with boundary 
conditions on Pi at 2; = and on P3 at z = d. (Pi and 
P3 are the densities of downward and upward directed 
chains, respectively.) From Eq. (|^), the stress components 
are simply 



4« 



= 2P2 + -(Pi 
= 0. 



(14) 

(15) 
(16) 



Before proceeding to the full solution of Eqs. (|3|), it is 
instructive to consider the linearized theory in the vicinity 
of P = 0: 



5. Pi 



n/3 



G; d,P. 



V3 



G. 



(17) 



Since G = Pi +P3 is a positive constant, one sees immedi- 
ately that P3 will become negative for sufhciently large z. 
In fact, for arbitrarily large values of P3(0), one sees im- 
mediately that P'i(d) must become negative for d > \/3/2, 
since the smallest possible value of G is P3(0). Thus the 
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linear theory produces unphysical results for systems of 
dimensionful depth A\/3/2 or greater. The problem can 
be traced to the divergence of chain densities mentioned 
above that occurs when the system is sufficiently deep 
to allow an appreciable number of force chains to "turn 
around" . 

In marked contrast to the linear theory, the full the- 
ory expressed in Eq. ([l^ ) admits physical solutions for all 
possible specifications of Pi{0) and Psid) for arbitrarily 
large d. Let B denote Pi{0)P3{d). The solutions are 



Pi 



Y — 7+ tanh 



S f - 7+ coth 



V3 



2 G7+ 
\/3 (1+G) 



f +7- 



tan 



G7- 



V3 (l+G) 



C) 



G > 2, B > 1 



G > 2, B < 1 



G < 2 



(18) 

where 7± = yJ±{G'^ ~ 4) /2 and G and C are constants 
determined by the specification of Pi (0) and P3 [d) . 

We will see below that a solution with non-negative 
densities exists for arbitrary choices of i-i(O) and P-i{d) 
and arbitrary values of d. On the other hand, for choices 
of P3(0) and/or Pi{d) within some ranges, the equations 
lead to negative densities or no solutions at all. This is a 
mathematical expression of the physically plausible state- 
ment that one cannot specify in advance the densities of 
chains directed outward from boundary. 

Exactly which specification of i-i(O) and P3{d) corre- 
spond to a given physical situation is not entirely clear. 
For a slab squeezed between two plates, however, a nat- 
ural assumption is that the boundary conditions should 
be symmetric: P^id) = Pi(0). In this case, C = d/2. 
Two particular solutions for symmetric boundary condi- 
tions are shown in Fig. ^ Part (a) shows the case G > 2, 
which yields extended flat regions in the upper and lower 
halves of the slab that are connected by a transition re- 
gion. The flat regions correspond to a uniform stress that 
is anisotropic: axx 7^ (^zz- The width of the transition re- 
gion is of order A, with an exponentially fast approach 
to the plateaus on either side. Surprisingly, the transition 
region is marked by a strong variation in the horizontal 
stress. Though the effect of such a transition region on 
the dynamical properties of the system is far from clear, 
it is interesting to note that the theory gives rise to an 
localized inhomogeneity that could influence shear band 
formation. The increased horizontal stress in the transi- 
tion region can be intuitively understood by considering a 
case of very large Pi(0) and P^id). In such a case, nearly 
all the chains in the top half of the sample are down- 
ward pointing, and nearly all the chains in the bottom 
half are upward pointing. The density of horizontal chains 
is determined by the splitting length, A. In the transition 
region, however, there is a high density of intersections of 
upward and downward pointing chains, which result in fu- 
sions that generate horizontal chains. Hence the density of 
horizontal chains is higher in the transition region, being 
determined by the distances chains propagate before fus- 
ing, rather than the longer distance required for splitting. 



5 
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Fig. 4. Typical solutions of the discrete Boltzmann equation 
with symmetric boundary conditions, (a) A generic solution 
with an exponentially localized transition region in the middle 
of the slab, (b) A solution with power-law boundary layers and 
a plateau in the middle, corresponding to Pi 4- -P3 = 2 ~ e. 



Fig. ^(b) shows density profiles for G < 2. In this case, it 
is the edges that deviate from the fixed point plateau in 
the middle, but the approach to the fixed point is a power 
law, ~ 1/z, so one cannot clearly define a boundary layer. 

The complete structure of the solution space is dis- 
played in Fig. ^. All points with P1P3 — 1, shown as a 
thick dotted line, are fixed points. Since Pi + P3 is a con- 
stant, each trajectory is a line with slope —1. The arrows 
indicate the direction in which the systems moves with 
increasing z. For trajectories in the region above the fixed 
line. Pi decreases with depth and P3 increases, and vice- 
versa for trajectories below (or to the left of) the fixed 
line. For small systems, trajectories in the lower left cor- 
ner, near the origin, may be relevant. For deep systems, 
however, the relevant trajectories must be ones that pass 
very close to some fixed point where Pi and P3 can stay 
nearly stationary for a long "time" . 

Let us now consider the possibilities for satisfying var- 
ious boundary conditions. Let U be the "starting" point 
(Pi(0), P3(0)) and V be the "ending" point (Pi(d), Psid)). 
First, we show that for any d, any Pi(0), and any Psid), 
U can be chosen in the first quadrant in such a way that 
V lies in the first quadrant; i.e., there exists a trajectory 
with no negative chain densities that satisfies the bound- 
ary conditions. To see this, first note that Pi{d) can never 
be negative. If U is chosen below the fixed line. Pi al- 
ways increases with increasing z. On the other hand, since 
trajectories can never cross the fixed line, if U is above 
the fixed line, the entire trajectory is confined to the first 
quadrant. To complete the proof, it is sufficient to show 
that as U varies in the first quadrant along the vertical line 
corresponding to any given choice of Pi(0), P:i{d) takes 
on all possible positive values. Consider U = {x,y). For 
any nonzero value of d, P3(ci) is negative for y = 0, since 
all trajectories beginning on the x-axis pass immediately 
into the fourth quadrant. As y is increased toward the 
fixed point y = Ps^d) rises continuously to l/x. As y 
is increased beyond l/x, Psid) increases without bound. 
Hence, given Pi(0) — x one can always choose y such that 
Psid) has any specified positive value. Fig. ||(b) indicates 
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with the assertion that physical constraints permit the 
specification of Pi{0) and /3(d) only. 

For deep systems (large d), there are three types of 
typical solutions, corresponding to the three expressions 
in Eq. (|TF 



P,(0) small, P/d) large _ 
\\\ ^Pj(O) small, Pj(d) small 




Fig. 5. (a) Trajectories of Eq. (|T|) in the P1-P3 plane. The 
curved line, P1P3 — 1 is a set of fixed points. The diagonal 
lines indicate the trajectories, with the arrows showing the 
direction of fiow with increasing z. For systems deeper than 
A, the typical length a force chain propagates before splitting, 
trajectories must either start very close to the line of fixed 
points (e.g. at point B), or on a trajectory that will take the 
system close to a fixed point (e.g. at point A), (b) The shaded 
regions indicate where trajectories must start in order to satisfy 
various possible boundary conditions for deep systems. The 
starting points must lie either in the vicinity of the fixed line, 
or just below the line of slope —1 tangent to the fixed line. 
For Pi (0) small and P3 (d) large, the starting point may be on 
either side of the fixed line. 



roughly where U must be chosen in order to satisfy vari- 
ous boundary conditions when d is large. For d small or of 
the order of unity, the trajectories generally lie completely 
within the unshaded regions of the diagram. 

Suppose now that we attempt to specify Pi and P3 
at the top boundary. In other words, we specify U com- 
pletely. In this case, problems arise for sufhciently large d 
whenever U is below the fixed line for Pi(0) > 1 or be- 
low the line Pi + P3 = 2 for Pi(0) < 1. Each trajectory 
in these regions passes into the fourth quadrant at some 
finite value of z, generating the unphysical result that P3 
becomes negative. Similar difficulties arise if one attempts 
to specify V rather than U. Finally, if one attempts to 
specify Pi{d) and ^3(0), there are no solutions at all for d 
large and Pi{d) > P3(0) + 1/-P3(0). The fact that certain 
regimes are prohibited in each of these cases is consistent 



Type I: Pi(0) > 1 and Pi (0)^3 (d) > 1. The motion must 
take place in the region above the line of fixed points, 
starting close to the line, at a point such as A in Fig. ||. 
For symmetric boundary conditions, the motion will 
also end near a fixed point, producing the behavior 
shown in Fig. ^(a). 

Type II: Fi(0) > 1 and Pi (0)^3 (d) < 1. The motion takes 
place in the lower right portion of the plane. The tra- 
jectory must start just below a fixed point. Such solu- 
tions occur if the bottom boundary condition is taken 
to be P3{d) — 0. Following such a trajectory backwards 
(up from the bottom of the sample) one finds the den- 
sities approach their fixed point values like 1/z. We 
will not pursue this case further here, but note that 
related asymmetric solutions might be relevant when 
gravity is included in the model. 

Type III: Pi(0) < 1 and Psid) < 1. For large d, the mo- 
tion must take place on a trajectory just below the 
one marked B in Fig. ||, so as to pass very close to 
the fixed point at Pi = P3 = 1. The density profiles 
associated with these trajectories are similar to those 
of Type II, but show increasing deviations from the 
fixed point values at both the top and bottom of the 
slab. Density profiles for a Type III solution are shown 
in Fig. ^(b). (Note the difference in horizontal scale 
between parts (a) and (b).) Though these solutions 
are mathematically consistent, they have the counter- 
intuitive physical property that the density of upward- 
pointing chains at the top surface is larger than the 
imposed density of downward-pointing chains. 

All other solutions are related to these three types by re- 
fiection through z — d/2 and interchange of Pi and P3. 

For deep systems, the system is primarily composed of 
regions in which (PijPs) lies close to a fixed point. The 
nature of the directed force chain networks at the dif- 
ferent fixed points is therefore of interest. First, we note 
that the networks associated with generic fixed points are 
anisotropic, having different densities for chains in differ- 
ent directions. Letting p designate the value of Pi at a 
fixed point, the value of P3 at that fixed point is 1/p. We 
emphasize that this is not a result of intrinsic anisotropy 
in the material, but rather a result of the boundary con- 
ditions. The situation is similar to that of an elastic ma- 
terial that may be isotropic when unloaded, but becomes 
anisotropic when subjected to a uniaxial stress. There is 
a crucial difference, however, in that the force chain net- 
work has no stable unloaded state. Thus, except under 
conditions of pure hydrostatic pressure, the response of 
the force chain network to perturbations in its loading 
will generally be anisotropic. The exceptional case is the 
fixed point at p = 1, where the system is indeed isotropic. 

Second, we may compute two distinct quantities that 
might be thought of as analogous to the Poisson ratio of 
standard elasticity theory. The first is the ratio Aa^x/ Aazz 
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obtained upon increasing the strengths of the forces in aU 
of the force chains. Since this does not change the network 
at aU, the resulting ratio is just equal to Oxxl <^zz at the 
fixed point. Using Eq. (|lj), we find 



(19) 



which varies monotonically from 1 at p = 1 to 1/3 for 
p — » 00. 

The other quantity analogous to the Poisson ratio is 
obtained by considering small increases in the applied den- 
sity of force chains, rather than their strengths. In this 
case, the system moves along the line of fixed points and 
we calculate 



1^2 



dpazz 



1/3. 



(20) 



It is amusing to note that a Poisson ratio of 1/3 arises 
also in studies of ball-and-spring networks with energies 
expressible as a sum of two-body central forces. We 
shall see in Section |]^, however, that this feature is spe- 
cific to the horizontal orientation of the 6-fold star vectors. 
For the vertical orientation, 1/2 varies with p. 

It is also interesting to note that Vi and V2 lie within 
the range of stability for 2D isotropic elasticity, in contrast 
to the Poisson ratio computed in Ref. jl^ from the linear 
theory of DFCN's. This is an indication that there are 
significant differences between the elasticity theories in the 
vicinity of a nontrivial fixed point and that obtained in the 
vicinity of the origin. 



4 Response to a localized force 

The response of the directed force chain network to a lo- 
calized force applied at the top boundary may be com- 
puted via linearization about the appropriate fixed point. 
In general, for a DFCN with a continuum of force chain 
intensities, there are two distinct ways to apply a localized 
perturbation: (1) by changing the strength of some of the 
force chains injected at the boundary; or (2) by adding 
some new force chains of chosen intensities. In the case 
where splitting and fusion angles are always 120°, how- 
ever, option (1) is not available. Because V's and V'/ permit 
interactions only among force chains of the same intensity, 
increasing the strength of one chain would be equivalent 
to removing that chain from the existing network and cre- 
ating an entirely new network completely decoupled from 
the original, which would lead right back to the problem 
of the failure of the linear theory in the vicinity of the ori- 
gin. Thus we must consider option (2), in which we inject 
a low density of new force chains with intensity /* along 
one or more of the directions . 

In the fixed point DFCN, the chain densities are trans- 
lationally invariant. This permits the derivation of decou- 
pled equations for the different Fourier modes of a pertur- 
bation applied at the top surface. We begin from the six 
equations obtained from Eq. (0): 



V3 



+ (PoP2-PrP3- AP5), (21) 
dxP2 = -P2 + P3 + P1 

+ (P1P3 - P2Pi - P2P0) , (22) 

— /3 1 

-^a.Pg + ^dxPs = -P3 + P2 + P4 

+ (P2P4 - - P3P1) , (23) 

and similar equations for i-4, P5, and Pq. Fixed point solu- 
tions of these equations have the form (Pf, P?, P^, Pq) — 
(p, 1/p, l/p,p), with P2" = P5O given by Eq. (|l|). Let p„ = 
P„(x, z) — P,j(a;, z) be the deviations from the fixed point 
and define Fourier coefficients (3„(g, z) via 



dq Qn{q,z) e" 



Linearization of Eq. (Ell) in Q yields 



dz 



Q(g) = M(q).Q(q), 



(24) 



(25) 



where Q is the column vector {Qi,Q3,Qi,Qo) and the 
complex elements of M are complicated algebraic func- 
tions of p and q. Analytic expressions for the eigenvalues 
of M may be obtained, but they are exceedingly compli- 
cated functions of p and q and we do not reproduce them 
here. (They involve solutions of quartic equations whose 
coefficients are polynomials of seventh order in p and third 
order in q.) The eigenvalues and eigenvectors of M(g) and 
M(~q) are related by the condition that the response to 
a real perturbation must be real. 

To compute the response to a localized perturbation 
applied at the top of a semi-infinite system, we seek solu- 
tions for which all QnS vanish as z — > 00. The symmetry 
of the Boltzmann equation and the fixed point solution 
under x —x guarantees that for every complex eigen- 
value K of M((7), there will be a complex conjugate partner 
K*. In addition, though the z — > — z symmetry is broken 
by the fixed point solution (since P^ ^ P3), there are 
two eigenvalues with negative real part for each q and two 
with positive real part. In constructing the perturbation, 
the coefficients of all eigenvectors having eigenvalues with 
positive real parts must be zero, since for these modes 
the Qn's grow exponentially with increasing z. By tak- 
ing appropriate linear combinations of the eigenvectors 
associated with the eigenvalues that have negative real 
parts, one can arrange to satisfy any specified boundary 
condition on pi{x^ 0) and pq{x, 0). Note that ^3(0;, 0) and 
P4(x,0) are then determined, which is consistent with the 
general rule that we may specify the densities of inward- 
directed chains only. 

For simplicity, we restrict the presentation here to ap- 
plied perturbations that are symmetric under refiection 
about a vertical line: perturbations with ^0(2;) — Pi{—x). 
The general form of the solution for p — (pi,P3,P4,Po) is 
then 



p(a;,z) 



dzPi + IdxPi = -Pi + P2 + Po 



o 2 
3 = 1 



Rc[Kj]js 



cos(ga: — liR[Kj]z) ej 



(26) 
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Fig. 6. The spectrum of eigenvalues of the linear operator 
governing perturbations in the vicinity of a fixed point for the 
case of the horizontal orientation of 6-fold star vectors, (a) The 
case Pi = 3. Thick (thin) lines are the real (imaginary) parts 
of K. A, B, and C mark regimes corresponding to different 
behavior of the response function. (See text.) (b) The case 
P\ = 1.001. Only the real part of k is shown. The quadratic, 
linear, and flat behaviors in the regimes A, B, and C are clearly 
visible. 



where Kj is an eigenvalue with negative real part, Gj is its 
associated eigenvector, and aj{q) is a constant determined 
by the boundary conditions. 

Fig. |^(a) shows the eigenvalues K{q) of M, computed for 
the case p — 3. The thick lines indicate the real parts of k 
and the thin lines the imaginary parts. The lower thick line 
and the two thin lines that emerge from the origin above 
and below the axis correspond to a complex conjugate pair 
of roots, the two roots with negative parts and therefore 
the only ones relevant for the response in a semi-infinite 
system. The spectrum may be divided, roughly, into three 
regions. In region A, corresponding to the smallest values 
of both q and |Re[K]|, neglecting all terms of order q"^ or 
higher in M we find 



K « —Dq + icq, 

where c and D are real constants: 

4 1 



(27) 



D = 



(28) 



V3 - p-^ ' 
1 

There is a transition region B with intermediate q and 
|Re[K]|, and finally, region C, a plateau in |Re[K]| for large 

q- 

The response at large z is dominated by the slowest 
decaying modes, those with |Re[K]| closest to zero. From 
the observation that there are no values of q other than 
zero for which Re[K] vanishes, it is clear that these are 
the l ong wavelength modes with the disperion relation of 
Eq. (27). Thus for large z we can approximate the integral 
in Eq. (pi) by 



dq e 



-Dzq^ 



-1/2 



[cos{q{x — cz))ei 

+ cos{q{x — cz))e2] (29) 
"ei-l-e 62 . (30) 



The response at large depths consists of two Gaussian 
peaks that propagate away from x = at an angle of 
tan~^ c and have widths proportional to \J Dz. 

Eq. (^) indicates that c corresponds to angles of prop- 
agation along the directions Vi and Vq for all of the fixed 
points. This appears to be special to the hor izontal ori- 
entation of the star vectors. (See Section 5.1 below.) D 
diverges as the isotropic fixed point, p = 1, is approached. 
This means that for nearly isotropic fixed points, the peaks 
at large z become increasingly broad. In addition, the re- 
gion over which the spectrum is quadratic (region A), be- 
comes smaller and smaller, so that the emergence of the 
two peaks occurs only for exceedingly large z. Fig. ^(b) 
shows the real part of the spectrum for the case p = 1.001, 
plotted on a logarithmic scale. In this case we see that the 
transition region B is a clearly defined regime of linear dis- 
persion. Linear dispersion generally implies either a single, 
centered peak in the response, or propagating peaks that 
broaden linearly with depth. 

The full response to a localized applied force can be 
computed numerically. Fig. ^ shows the response func- 
tion for the case p = 3. The curves were obtained by ap- 
proximating the integral over q in Eq. ( p6| ) with a sum 
over discrete values q — {n — l/2)Sq, with Sq — 0.01 
and 1 < n < 1000. This range of q's contains enough 
high g's to construct a relatively sharp initial Gaussian 
at z = 0, and enough low g's to observe the response at 
large z. For each q, the eigenvalues and eigenvectors of 
M are determined. We then obtain the linear combina- 
tion of the two stable eigenvectors required to produce 
the vector Q = {l,y,y,l), where y is unspecified. (The 
two conditions pi (g) = 1 and po (q) — 1 determine the lin- 
ear combination of the two eigenvectors completely, and 
hence determine ps = P4 = y, and p2 and in turn.) 
These linear combinations of eigenvectors are then multi- 
plied by a factor (oj = exp(— g^/10)) so that the sum in 
Eq. (26) at z = yields centered Gaussians for pi and 
which constitute the applied perturbation. The response 
is then determined at various values of z by summing the 
integrand of Eq. ( p6|) over the discrete set of q's. 

Fig. 0(a) shows the p„'s at the top surface. The iden- 
tical, centered Gaussian pi and po represent the imposed 
boundary condition. The other p„'s at the top surface are 
part of the response, as determined by the condition that 
P3 and p4 must vanish at z = oo. Part (b) shows profiles 
of along horizontal slices at various values of z. The 
curves have been scaled so that the trend toward diffu- 
sively broadening, linearly propagating peaks is evident. 
At large z, the peaks in this figure would become increas- 
ingly sharp while maintaining their same amplitude and 
position. Part (c) shows the relation of the various stress 
components at large z. 



5 Other choices of star vectors 

5.1 6-Fold star in the vertical orientation 

After the horizontal orientation of the 6-fold set of star 
vectors, the next simplest case is the vertical orientation 
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Fig. 7. Response functions for 6-fold star vectors in horizontal 
orientation, (a) Pn's at the top surface. Pi and Pq are identical 
and are specified as boundary conditions. The other P„'s are 
part of the response. Due to the symmetry of the boundary 
conditions, P4(x) = P3(— i) and Psix) = P2{—x). (b) Scaled 
profiles of azz for several values of z. From top to bottom, the 
curves correspond to 2 = 0.2,0.5,1,2,5, and 10. (c) Different 
components of the stress tensor at 2 = 10, corresponding to 
the thick line (b). Only positive x is shown. The symmetry of 
the boundary conditions implies Ozzix) — azz{—x), axx{x) = 
c^xxi-x), and axz{x) = ~Oxz{~x). 



of the 6-fold star. Calculations for the vertical orientation 
are an important means of testing the robustness of sev- 
eral features of the solutions found above. For horizontally 
uniform systems with reflection symmetry under x ^ —x, 
the Boltzmann equation reduces to a set of four coupled, 
nonlinear, ordinary differential equations for the variables 
Fo, A, P2, and P3. (See Fig. |.) Using the fact that azz 
is conserved, these can be reduced to a set of three. Exact 
analytical solution for the chain densities as a function of 
depth for generic boundary conditions is difhcult. Never- 
theless, it is easy to find a line of fixed points, which is 
all that is required for carrying out an analysis of the re- 
sponse as performed above for the case of the horizontal 
6-fold star. The fixed points identified below were found 




Fig. 8. A typical solution of the discrete Boltzmann equation 
with symmetric boundary conditions for the vertical orienta- 
tion of 6-fold star vectors. (Compare to Fig. ^). The boundary 
conditions are Po(0) = 4, Pi(0) = P5(0) = 2.5, Paid) = 4, 
P2(d) = P4,{d) — 2.5, with d = 10. The transition region in 
the middle creates a bump in axx- Boundary layers are also 
present. 



by making the assumption that the densities of chains in 
opposite directions are related by multiplicative inversion; 
i.e., P0P3 = P1P4 = P2P5 — 1- (This was inspired by the 
horizontal case, in which exactly analogous relations were 
found to hold, but we have not ruled out the possibility 
that other fixed points exist.) 

Numerical solution of the two-point boundary value 
problem yields solutions of the type shown in Fig. ^. The 
basic features of the solution are quite similar to those 
for the horizontal orientation of the star vectors. In par- 
ticular, there are flat regions corresponding to fixed point 
solutions and a transition region in the middle of the slab 
with enhanced horizontal stress. In the vertical case, how- 
ever, boundary layers at the top and bottom appear as 
well. These arise because the additional dimensions in P- 
space allow for trajectories that begin away from the fixed 
line, approach it rapidly, then eventually diverge from it 
again; i.e., the fixed line has both attracting and repelling 
directions. The precise forms of the chain density profiles 
in the boundary layers and in the transition region are 
sensitive to the details of the imposed boundary condi- 
tions, but the enhancement of a-xx in the transition region 
is a robust feature. 

The fixed points of Eq. (^) with x —x symmetry 
have 

p; P2 



Po 

which yield 



p'; Pi 



p 



P:^ 



P 



3(p + p-i) 



3 



l + 4(]3 + p-i) 



(31) 

(32) 
(33) 



for the two types of "Poisson ratio" . Here, as for the case 
of the horizontal orientation of the 6-fold star vectors, vi 
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Fig. 9. The spectrum of eigenvalues of the linear operator 
governing perturbations in the vicinity of a fixed point for the 
case of the vertical orientation of 6-fold star vectors. Only the 
eigenvalues with negative real parts are shown, (a) The case 
Pi = 2. Thick (thin) lines are the real (imaginary) parts of k. 
(b) The case Pi = 1.1. In both (a) and (b), the magnitudes of 
the imaginary parts have been multiplied by 15 for clarity. 



varies monotonically from unity to as p ranges from 
1 (the isotropic fixed point) to oo. Unlike the horizontal 
case, however, V2 is not constant. It is again 1/3 at the 
isotropic fixed point, but decreases monotonicaUy to zero 
for large p. 

To determine the response function, we work in the 
vicinity of a generic fixed point, neglecting the effects pro- 
duced by the boundary layer. Fourier analysis of the lin- 
earized equations yields the spectrum shown in Fig. ^. The 
spectrum at low q has the form of Eq. ( p7[ ) , with 

^ ^ 6(2r4+ + n+ + 8r2+ + 9ri+ + 14) 
(4ri+ + if (2r3_ -I- r2_ - 2ri_) 
c = V3(4ri+ + l)-i. (34) 

where we have defined rn± = ± For strongly 

anisotropic fixed points (i.e., large p) we have c ~ (3/8)p^^ 
and D « (4p/3)~^/^, implying, at large z, slowly dif- 
fusing peaks propagating close to the vertical direction. 
For nearly isotropic fixed points p = 1 + e, we have c ~ 
y/lji (1 - 2eV9) and D « (l/3)e-\ implying rapidly 
broadening peaks propagating along directions close to 
halfway between the star vectors. Thus, in constrast to 
the case of the horizontal star orientation, the direction 
of propagation of the peaks varies continuously from 30° 
near the isotropic fixed point to near 0° at the strongly 
anisotropic fixed points: the direction of propagation of 
the peaks is not tied to the star vector directions. 

As in the case of the horizontal orientation, the region 
in q space over which the dispersion relation is quadratic 
shrinks to zero as the isotropic fixed point is approached. 
We also note that near the isotropic fixed point, where D 
becomes large, we can eliminate e to obtain c ~ -^/l/S (1 — 
(2/81)Z)~^), whereas for the strongly anisotropic fixed 
points, where D is small, we can eliminate p to obtain 
c « 1)^/2. These relations between propagation direction 
and diffusion rate differ markedly from relations obtain 
from models in which the characteristic propagation direc- 
tion and diffusion rates are both controlled by a parameter 
corresponding to the strength of the disorder, where one 




-1 -0.5 0.5 1 -1 -0.5 0.5 1 




X X 



Fig. 10. Response functions for 6-fold star vectors in verti- 
cal orientation. Solid lines indicate CTzz and dashed lines a^x- 
Thick lines correspond to a boundary condition in which Pq is 
a Gaussian of unit strength and Pi and P5 are 0. Thin lines 
correspond to a boundary condition in which Po = and Pi 
and P5 are Gaussian of strength 2, so as to make (Jzz{z — 0) 
the same for both cases. Note the difference in both horizontal 
and vertical scales in the different panels. 



finds c « Co - aD" with a > 0. [|jT^ This suggests that 
the DFCN is not in the weak disorder regime. 

Typical response functions are shown in Fig. |l^. Two 
independent choices for the application of a localized ver- 
tical load are shown: 

( (e,0,0)exp(-(a;/0.09)2) 
{Po,Pi,P5) = < and . (35) 

[ (0,2e,2e) exp (-(a;/0.09)2) 

These two choices produce the same profile for azz{0), 
but different profiles for axx{0), as shown in (a). At the 
small depth shown in (b), the two give different responses 
with somewhat complicated structures. At the interme- 
diate depth shown in (c), the responses differ as well, 
but both have something like the Lorentzian structure fa- 
miliar from the response of semi-infinite, standard elas- 
tic materials. pO| In particular, both show a central peak 
with long tails. It is not until we reach depths signifi- 
cantly larger than A (which has been scaled to unity) , as in 
(d), that we see the two diffusively broadening Gaussian 
emerge and the responses become quite similar. 

In the asymptotic regime at large depths, Fig. |l^ sug- 
gests that a XX becomes proportional to azz- In fact, we 
can show both that axx ~ c^cFzz and that c — ^JT^ by 
considering the eigenvectors at low q associated with the 
dominant branch of the spectrum. Let 

axM = S^'^ •e(g), 

a,4q) = C(2).e(g); (36) 
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where C'^' and S^^^ are row vectors with components 
cos^(7rj73) and sin^(7rj73), respectively, and e{q) is the 
eigenvector of M{q) associated with the dominant branch. 
e(g) converges as q approaches zero precisely to the 9 = 
mode corresponding to shifts along the fixed line. But 
these shifts are just the ones used to compute 1^2 ■ Thus 
as z increases and the surviving modes have q^s closer to 
zero, and the ratio (Jxx{m) / f^zzio) for all of those modes 
approaches V2. To see that the ratio must equal c^, note 
that on large length scales the response appears as two 
narrow lines that propagate along the directions x = ±cz. 
Newton's laws require that the stress along these lines 
correspond to forces directed along the line, which imme- 
diately implies <7xx = c^f^zz along the line. 

The crossover from a roughly standard elastic form to 
a hyperbolic one with diffusive broadening is an unusual 
feature. We note that linear elasticity theories, even in 
the unstable regime where the equations are hyperbolic, 
always yield scale invariant response functions for a semi- 
infinite medium. The spectrum is always purely linear, so 
that the response is a function of x/ z only. In the DFCN, 
however, the spectrum has a richer structure and the re- 
sponse at large z may appear quite different from that 
at small or modest depths. Though the discrete model 
studied above does not represent a realistic picture of a 
disordered granular material, the result does suggest that 
experimental data should be interpreted cautiously: the 
appearance of a response with a typical elastic form may 
be a deceiving effect of working with a relatively small 
sample. 

Fig. 11 shows a contour plot of the response for the sec- 
ond choice of perturbation in Eq. (^5|). The figure clearly 
shows the emergence of the double-peaked structure at 
large z from a single peak at intermediate z. It also shows 
multiple peaks near z = that have not yet been men- 
tioned. These peaks, which are of very low amplitude, form 
and decay rapidly. The feature in the spectrum that is re- 
sponsible for them is the broad peak in the uppermost 
branch of the real part in the neighborhood oi q — 2.5. 
(See Fig. |.) Modes near this value of q decay slightly 
more slowly than others nearby. 



5.2 n-fold symmetric stars in 2D 

For sets of star vectors with the symmetry of a regular 
n-gon, as depicted in Fig. |l^, the full nonlinear problem is 
much more difficult due to the fact that splitting generates 
chains with different intensities. Nevertheless, analysis of 
the Boltzmann equation with fusion terms neglected can 
be carried out and yields useful insights. 

When the fusion terms in Eq. (H) are neglected, the 
equation can be multiplied by / and integrated over / 
to obtain a new Boltzmann equation governing the spa- 
tial variations in the quantity F{9) = J^df fP{f,6), the 
force density associated with chains in the 9 direction. | p6[ 
This equation can be analyzed in detail for the case where 
all chains are directed along the discrete set dj = 2TTj/n 
for < j < n, depicted in Fig. |l^. Consistency with the 




Fig. 11. Contour plot of a^z for the vertical orientation 
of 6- fold star vectors an applied chain density (po,Pi,P5) = 
(0, 2e, 2e) exp( — (s/O.OQ)^). The contour hnes are at the set of 
values (-0.001, 0., 0.001, 0.005, 0.007, 0.01, 0.014, 0.02, 0.028, 
0.04, 0.056, 0.08). The emergence of the double-peaked struc- 
ture from a single peak between z = 0.5 and 4.0 is clearly 
visible. For z < 0.5, a number of small, rapidly decaying peaks 
appear corresponding to wavenumbers near the broad peak in 
the spectrum near q = 2.5, which is visible in Fig. H. 




Fig. 12. The n-fold star used for analysis of the linear theory. 
Each angle is 27r/n, where n is an integer not divisible by 4. 



restriction to this discrete set can be maintained by as- 
suming that splittings are always symmetric and outgoing 
chains make angles of 27r/n with the incoming one. (For 
n > 8 one can also have splittings at angles of An/n, and 
so on, but the simplest case is sufficient for our purposes.) 
To avoid the need for special treatment of the horizontal 
directions, we also consider only n's that are not multiples 
of 4. 

Let Fj be the force density for chains in the direction 
9 = 27rj/n. After rescaling of A to unity, the continuum 
equations for the F's are 



cos(j0) d,Fj + sin{j9)dxFj = -F, + — (F,_i + ^;,-+i) , 

(37) 

where 9 = 27r/n, ci = cos0, and indices are taken modulo 
n. An analysis of this equation in the slab geometry is pre- 
sented in the appendix. It shows that the lack of solutions 
to the linear problem in the 6-fold case is generic: force 
densities always diverge when the system is sufficiently 
deep, even though the intensities of the chains become 
exponentially small after many splittings. 
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5.3 Three dimensions 

In 3D, the analogue of the 120° vertex for splitting and 
fusing is the tetrahedral vertex. Three types of events are 
possible: splitting, in which one chain splits into three; fu- 
sion, in which three chains fuse into one; and scattering, 
in which incoming chains along two of the tetrahedral di- 
rections produce two outgoing chains along the remaining 
two directions. In principle, scattering events should be 
more common than fusions, which require terms in the 
Boltzmann equation of order P^. Scattering events do not 
change the total force density in the system, however, and 
therefore are not sufficient to regulate the divergence in 
the linear theory. Though scattering events should clearly 
be included in the Boltzmann equation for completeness, 
fusions are the essential process. 

A suitable star of vectors is formed from the 20 vec- 
tors pointing to the faces of an icosahedron. (The appar- 
ently simpler choice of the star of vectors pointing to the 
faces of an octahedron turns out to possess certain sym- 
metries that lead to nongeneric behavior.) Each vector on 
the icosahedral star is a part of two different tetrahedra. 
One can assume, for example, that the splitting of a chain 
occurs along each tetrahedral set with probability 1/2. 

Choosing the orientation of the icosahedral star to in- 
clude one vector in the downward vertical direction, and 
neglecting scatterings, at least one line of fixed points 
can be found. In this orientation, there is one vector v 
with z-component Sz = 1, three with = \/5/3, six 
with Sz = 1/3, and 10 related to these by inversion. The 
fixed points can be parametrized by p, with the densities 
of chains as follows: P = p^/\/2 for all — 1 chains; 
P = pV\/2 for aU Sz = V^/B chains; P = p/V2 for aU 

= 1/3 chains; and for all chains related to these by in- 
version, P 1/(2P). This gives for the "Poisson ratios" 



V2 



8ri+ + 27-2+ 



2ri+ + 10r2+ 3r3+ 

8ri_ -H 4r2- 
2ri_ + 20r2- + Org- 



(38) 
(39) 



v\ varies from 2/3 at the fixed point, which is outside 
the range of stability for a standard, isotropic elasticity 
theory, to at large p. vi varies from 16/69 to and thus 
is always within the stable regime. 

This preliminary treatment of a 3D model suggests 
that there is nothing fundamentally different from the 2D 
models studied above. The spectrum and response func- 
tions in the vicinity of a fixed point can be calculated nu- 
merically. Extracting the asymptotic form of the response, 
however, is not as straightforward as in the 2D case. Thor- 
ough analyses of the nature of the response and the effects 
of the scattering terms are beyond the scope of this work. 



5.4 Long-wavelength theories for a continuum of chain 
directions 



to a discrete set of directions. It was shown that for the 
linear Boltzmann equation, multiplication by / and inte- 
gration over / could be carried out to obtain equations for 
the force densities F(ff) = df fP{f,9). These in turn 
could be cast in terms of equations for the stress tensor 
under the assumption that F{9) was nearly isotropic. This 
assumption is clearly not true on length scales of order A, 
but was argued to be true at large length scales. In terms 
of F, the following three quantities were defined: 



1 

P = n 



dQF{e)- 



3 = -J dQF{9)n- 
I d[2F{9) n(g)n, 

J —TT 



(40) 



where df2 = d9 /2tt as defined above. Here p is the hydro- 
static pressure, a is the stress tensor (see Eq. (||), and J is 
a vector field that plays a similar mathematical role to the 
displacement field in standard elasticity theory. Under the 
assumption that F{9) is nearly isotropic, one can express 
F(n) as a linear combintation of fi • J, and n • cr • fi, 
then use Newton's laws and a compatibility condition on 
J to derive a partial differential equation for a. That equa- 
tion turned out to be elliptic, which suggested a response 
similar to that of an elastic material, 

From the results of the present work, it appears that 
the assumption of the isotropy of F{9) at large length 
scales was not justified. Analysis of the linear theory, as 
discussed in the appendix, shows that the isotropic fixed 
point F = is unstable to perturbations that favor large 
F{9) for 9 near 7r/2, and that the nontrivial fixed points 
are generically anisotropic. Thus the analysis presented 
m Ref. @ should be generalized to the anisotropic case. 
It is not clear whether such a generalization can preserve 
the close correspondence between J and the displacement 
field of standard elasticity theory. 

In addition to the issue of anisotropy, this program en- 
counters difficulties related to the neglect of the nonlinear 
terms in the derivation of the equations for F{9). Naive 
attempts to carry out the same program without neglect- 
ing the nonlinear terms fail: the structure of the (j)f terms 
in Eq. (||) make it important to calculate the full distri- 
bution P{f,9). It is not possible, without making further 
approximations, to recast the equations in terms of F(9) 
alone, m 

It should also be noted that the implicit passage to 
the limit of large length scales corresponds to keeping only 
terms of order q in the expressions for c and D. Such a cal- 
culation does not pick up the quadratic dependence of D 
on q. Thus, in these theories, whenever the response is hy- 
perbolic the apparent value of D is zero, correponding to 
delta-functions that propagate without broadening. Inclu- 
sion of the terms in the expansion of F containing nonzero 
wavevectors q would be required in order to compute D. 



In Ref. ||l6|, equations for the stress tensor were developed g DisCUSSion 

for a directed force chain network without the restriction 
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6.1 Numerical simulations 



6.3 Boundary conditions 



There is no known algorithm for numerically generating 
configurations that are consistent with the Boltzmann equa- 
tion. Simulations previously reported were done on 
small systems and used a regularization scheme that is 
uncontrolled. The difficulty encountered in these simula- 
tions is that the effects of fusions can create "causality" 
problems: a chain that gives rise to many others through 
multiple splittings may itself be altered by a fusion with 
one of its offspring. For small systems, this problem does 
not occur and the pseudo-elasticity theory derived from 
the linear Boltzmann equation appears to give a good de- 
scription response.]!^ For large systems, however, it is 
an important part of the physics. The algorithms used in 
Ref. [|l6| for generating DFCN's grind to a halt when the 
system size is increased to a few times the length A, due 
to the generation of infinite loops associated with fusions. 

One approach to avoiding the causality problem would 
be to neglect fusions altogether. To avoid an exponential 
explosion in the number of chains (and the total force den- 
sity), one must then introduce a lower cutoff on force chain 
intensities. All chains with intensity lower than the cutoff 
are simply ignored. Again, this approach may be reason- 
able for small systems where the chains being neglected 
are just those that were generated in splitting events where 
one daughter chain was in nearly the same direction as the 
parent and the other was very weak. For very large sys- 
tems, however, chains split many times before reaching the 
bottom. At large depths, the number of chains surviving 
above the cutoff will eventually decay to zero, so no stress 
at all will be transmitted to large distances. The appar- 
ently justifiable neglect of small forces leads to substantial 
violations of Newton's laws. 

It appears that the numerical generation DFCN's will 
require algorithms for relaxing candidate configurations in 
a nontrivial way so as to arrive at networks with statistical 
properties consistent with the Boltzmann equation. 



6.2 Strong vs. weak disorder and the importance of 
discrete directions 

Hyperbolic systems with weak disorder yield diffusively 
broadening, propagating peaks. Q The strong disorder re- 
gime has not previously been accessible. In one respect, 
our discrete model exhibits strong disorder: splitting events 
cause chains to rapidly lose memory of the direction of 
their ancestors. In another sense, though, the disorder may 
be weak: the orientation of the star vectors exhibits no 
fluctuations. In any case, the discrete model may apply 
rather directly to the case of a hexagonal array of disks 
with a random distribution of vacancies. 

We are currently investigating the behavior of the Boltz 



A major unanswered question about the DFCN theory is 
how to determine the appropriate boundary conditions as- 
sociated with a particular physical situation. Assume, for 
example, that we clamp a slab of granular material be- 
tween to nominally flat pistons. How do we decide whether 
we are injecting a high density of weak chains or a low den- 
sity of strong chains? The problem is complicated by the 
fact that some part of the pressure on the pistons derives 
from the response to the injected chains, so our bound- 
ary condition must be chosen so as to produce the correct 
pressure at the associated fixed point; we cannot simply 
inject chains corresponding to the desired loading of the 
surfaces, though in the slab geometry we can always scale 
all force chain intensities by a common factor to produce 
any desired total vertical force on the top and bottom 
boundaries. 



6.4 Anisotropy and gravity 

As one might expect on symmetry grounds, generic fixed 
point solutions of the Boltzmann equation in the slab 
geometry are anisotropic even though the equation does 
not contain intrinsic symmetry breaking terms, and this 
anisotropy has a qualitative effect on the response. The 
source of the differences in response are traceable to the 
inherent nonlinearity of the system. The response near a 
given fixed point cannot be construed as the superposi- 
tion of the response in an isotropic system together with 
a homogeneous anisotropic deformation. 

Inclusion of gravity in the DFCN theory may be ac- 
complished by assuming that every grain is a source of a 
new vertical chain, which may or may not instantly fuse 
with other chains passing through the same grain. This 
would require the addition of a source term to the Boltz- 
mann equation that breaks the isotropy of the equation 
itself, not just the solutions with anisotropic boundary 
conditions. Similar considerations would apply for systems 
with an anisotropic fabric tensor; i.e., systems in which 
certain directions are favored by an intrinsic anisotropy in 
the packing geometry. Analysis of systems with intrinsic 
anisotropy is beyond the scope of this work. 



6.5 Standard elasticity 

Materials described by standard, linear elasticity theory, 
even if isotropic, always have Fourier spectra with expo- 
nents linear in q. Such theories can produce propagating 
peaks, as have been observed for ball-and-spring models 
with strong anisotropy p2[ . But the peaks must always 
broaden linearly, never diffusively. For classically stable 



mann equation for more general splitting and fusing functionsn |^] erials, the peaks will always broaden linearly with 



To make predictions for real granular materials that have 
not been carefully prepared to have orientational order, 
it is crucial to know whether the hyperbolic response at 
large depths survives in the absence of a globally specified 
set of favored directions. 



depth. For a semi-infinite slab, there is no parameter in 
the theory with dimensions of length, so the response is 
always a function of x/z only. Outside the domain of sta- 
bility, however, the response may be hyperbolic. We are 
currently investigating anisotropic models in this regime 
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in an effort to determine whether they can describe the 
bahavior of DFCN's at very large length scales. pi| 



6.6 Interpreting experiments 

The emergence of a two-peaked, hyperbolic response at 
large depths in the models we have studied is quite re- 
markable, particularly in the vertical orientation of the 
6-fold star, where the directions of propagation are not 
simply related to the star vectors. In a strongly scatter- 
ing system such as this one, where splittings cause large 
changes in the directions of force chains, one might naively 
expect any peaks to be strongly broadened at large depths, 
which would lead to a single centered peak in the response 
function. Indeed, this is what we observe for moderate 
depths: a system may have a two-peaked response close 
to the surface, but those peaks decay or broaden rapidly, 
and a single peak emerges. Surprisingly, in the AY-model 
two peaks that really are associated with hyperbolic re- 
sponse emerge at even greater depths, and the directions 
of propagation of these peaks are generally not the same 
as those of the original two peaks near the surface. 

Experiments on 2D systems have shown both propa- 
gating peaks and single centered peaks. 1 23, 2^,^ In in- 
terpreting these experiments it is important to note three 
things: (1) Correlations in the positions of grains at the 
surface with the layer just below may affect the interpre- 
tation of the distribution of chain directions injected at 
the surface. It would not be surprising to see structure 
in the response close to the surface, possibly even double 
peaks due to the immediate splitting of a vertical force 
chain. (2) For well-ordered systems, the splitting length 
A may be large enough so that the entire system is in 
the very shallow regime. Thus the two peaks observed are 
not necessarily indicative of a hyperbolic response at large 
scales. (3) Even for disordered systems, one may have to 
probe systems an order of magnitude or more larger than 
A before the true large system behavior becomes appar- 
ent. This may be especially difficult because the peaks 
decrease in amplitude at large depths, but the fact that 
diffusively broadening peaks become cleanly separated at 
large depths offers some hope that they could be resolved 
in future experiments. 
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A Failure of the linear theory near the origin 

In this appendix we present an analysis of the linear Boltz- 
mann equation for a discrete set of chain directions and 
symmetric splittings. The main purpose is to clarify the 
nature of the divergence and hence the failure of the linear 
theory. 
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Let Fj be the density of force chains with direction 9 = 
27rj/n, where n is an integer that is not a multiple of 4 
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and 6 = corresponds to the positive z direction. (See 
Fig. [l^.) Assume that sphttings always occur at an angle 
9 = 27r/n; i.e., when a chain in the j direction splits, the 
only option is to produce chains in the j — 1 and j + 1 
directions. Eq. (|37|), reproduced here, governs the force 
densities. 



cos{j9) d,Fj + sinijd) d^F, = -F, 



F. 



j+i) 



(41) 

We seek solutions to this equation in the slab geometry 
with a load that is uniform in the horizontal direction. 

Translational invariance ensures that the solution will 
be uniform in the horizontal direction, so the dx terms 
vanish. Letting Cj stand for cos{j9), we have 



d,F, 



In matrix form. 



1 

—F, 



1 



2cic. 



{Fj-i+F,+i). 



d,F=M- F, 



(42) 



(43) 



where F is the n-dimensional vector with elements Fj and 



M 



2cico 







2ciCo 



2cici 




2cic„-i 



2cj.ci 



2ciC2 C2 2ciC2 











— / 



(44) 

M is tridiagonal (except for the entries in the corners), 
but not symmetric. Its eigenvectors are not orthogonal, 
nor do they span the full space. One can check that both 



and 



C(2) • M = 



X(2) . M = 0, 



(45) 



(46) 



where C^^) and X^^) are the (row) vectors with compo- 
nents c| and Cj sm{j9), respectively. This implies that C*^^^ 
and X^^) are orthogonal to every eigenvector of M with 
a nonzero eigenvalue, and therefore that no matter how 
those eigenvectors grow or shrink with increasing z, the 
projection of their sum onto C^'^-' and X^-^^ will be con- 
stant determined by the amplitudes of the eigenvectors 
with eigenvalue zero. These projections corresponds to azz 
and axz, and their invariance is a direct consequence of 
Newton's laws. 

The matrix M can be block diagonalized to separate 
subspaces that are symmetric or antisymmetric with re- 
spect under j — > n — j. For present purposes, it is suf- 
ficient to consider only the symmetric subspace. Let 
denote the m x m matrix corresponding to the symmet- 
ric block. For n odd, we have m = (n -I- l)/2 and for 
n even m = (rt/2) + 1. has two eigenvalues equal to 
zero. One is associated with the eigenvector C with ele- 
ments Fj = Cj. The other is associated with a vector Ei 
that satisfies M • Ei = C rather than the usual eigenvalue 
equation M • E = AE. One can check that the solution 



is El = 5(1, 1, 1, • • •) — aC(^\ with a = ci/(ci — C2) and 
b = a(l — ci). 

The solution of the differential equations for an arbi- 
trary initial condition within the symmetric subspace is 



F = aoC + ai (Ei + zC) 



k=2 



E. 



(47) 



where the a^'s are determined by the boundary condi- 
tions. We emphasize the following three features of the 
solutions. First, projections onto C are not conserved as 
z varies, even though the eigenvalue associated with C 
is zero, because C appears also multiplied by z. Second, 
symmetry under z —^ —z and Fj — > Fn/2-j guarantees 
that if A is an eigenvalue of M^, —A is also, implying that, 
aside from the eigenvectors with zero eigenvalues, half of 
the eigenvectors grow exponentially with increasing z. Fi- 
nally, all the Efc's except Ei must have both positive and 
negative components, since their projection onto the pos- 
itive definite C^^^ vanishes. 



A. 2 Negative force densities and their origin 

To avoid negative force densities in a deep system, the 
boundary values of F must have zero projection on all 
eigenvectors with positive eigenvalues. Since azz = C^^^ • 
F, any boundary condition at z = that has nonvanish- 
ing azz must have oi ^ 0. But this implies that the zC 
term in Eq. (^7|) cannot be avoided. Since C has negative 
components, and all the contributions from eigenvectors 
with negative eigenvalues decrease with increasing z, the 
solution will always produce negative force densities at 
sufficiently large z. 

The appearance of negative force densities may be 
counter-intuitive. After all, the linear equations describe 
physically allowable processes that can never generate neg- 
ative densities of chains and never yield negative intensi- 
ties for any chain. The source of the negative -F's in the 
solution is a divergence that occurs when the system is 
sufhciently deep to allow enough splittings that an appre- 
ciable fraction of the descendants of a downward-pointing 
chain are directed upwards. This causes an exponential 
explosion in the number of chains in the system and a 
consequent divergence of the force densities. Even though 
the intensities of the descendants of a given chain decrease 
exponentially with generation number, the sum of their in- 
tensities diverges because in every splitting event from a 
chain with intensity /i to two with /2 and /a we have 
(/2 -I- fd,)/ fi = sec{9) > 1. The effect on the solutions to 
the differential equations is roughly analogous to the con- 
tinuation of J2m=o'''"^ ~ 1/(1 ~ beyond its radius of 
convergence r = 1, where it appears that an infinite sum 
of positive terms gives a negative result. In other words, 
the theory is not well defined for slabs of thickness larger 
than a "persistence length" for chain directions. 
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A. 3 The question of anisotropy at large scales 

The problem of negative densities in the hncar theory 
shows that the response of the DFCN system can only 
be properly defined for a pre-stressed system. This still 
leaves open the possibility that the same linear theory 
might describe deviations from some uniformly stressed 
state. 

Here we point out that in the slab geometry with uni- 
form horizontal loading, the instability of the trivial solu- 
tion, = 0, to the linear equations leads to anisotropic 
F^s for a broad class of sphtting functions ips- For discrete 
n-fold systems in which symmetric splittings (02 = ^s) are 
allowed with uniform probability for all angles up to = 
7r/3, numerical diagonalization of the matrix M of Eq. ( [44|) 
shows that the most strongly unstable q = modes cor- 
respond to eigenvectors E{9) peaked near 6 = ±7r/2. The 
eigenvalues grow more unstable and the eigenvectors more 
sharply peaked as n increases toward the continuum limit. 
Analytical arguments have also been obtained for certain 
special choices of V's(02, ^a)- 

The anisotropic nature of this instability casts doubt 
on the validity of the conjecture that the DFCN should 
appear isotropic on large scales. The isotropic DFCN does 
not appear to be an attractor for any but a perfectly tuned 
set of initial conditions. 



